Phase diagram of van der Waals-like phase separation in a driven granular gas 
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Equations of granular hydrostatics are used to compute the phase diagram of the recently dis- 
covered van der Waals-like phase separation in a driven granular gas. The model two-dimensional 
system consists of smooth hard disks in a rectangular box, colliding inelastically with each other 
and driven by a "thermal" wall at zero gravity. The spinodal line and the critical point of the phase 
separation are determined. Close to the critical point the spinodal and binodal (coexistence) lines 
are determined analytically. Effects of finite size of the confining box in the direction parallel to the 
thermal wall are investigated. These include suppression of the phase separation by heat conduction 
in the lateral direction and a change from supercritical to subcritical bifurcation. 
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I. INTRODUCTION 



Granular gases (gases of inelastically colliding macro- 
scopic particles) exhibit a plethora of symmetry-breaking 
instabilities and clustering phenomena . Their investi- 
gation is useful both for testing and improving the models 
of granular flow, and for a deeper understanding of far- 
from-equilibrium dynamics in general. In this work we 
focus our attention on a recently discovered phase sep- 
aration instability which occurs in a prototypical two- 
dimensional granular system: an assembly of monodis- 
perse hard disks in a box, colliding inelastically with each 
other and driven, at zero gravity, by a rapidly vibrating 
or thermal wall. An immediate consequence of the in- 
elasticity of the particle collisions is the formation of a 
laterally uniform cluster (the stripe state) at the wall 
opposite to the driving wall 0, 0- Both granular hy- 
drodynamics and direct molecular dynamic simulations 
show that this simple clustering state can exhibit spon- 
taneous symmetry-breaking instability leading to phase 
separation: coexistence of dense and dilute regions of the 
granulate (droplets and bubbles) along the wall opposite 
to the driving wall 0, H H 0, H H H3 . This far-from- 
equilibrium phase separation is strikingly similar to the 
gas-liquid transition in the classical van der Waals model. 
The objective of this work is a systematic investigation of 
the steady states in this system and computation of the 
phase diagram starting from the Navier-Stokes granular 
hydrodynamics. Granular hydrodynamics is expected to 
be an accurate leading order theory when the mean free 
path of the particles is much less than a characteristic 
inhomogeneity scale of the problem, and the mean time 
between two consecutive collisions of a particle is much 
less than any time scale the hydrodynamics attempts to 
describe. In addition, we should work with sufficiently 
low particle densities, when an account of binary colli- 
sions and volume exclusion effects is sufficient [Tl|. As 
will be shown below (see also Ref. Q), the requirement 
that the mean free path be small compared to the in- 
homogeneity scale implies in this system that particle 
collisions must be nearly elastic. 



For steady states with a zero mean flow, granular hy- 
drodynamic equations reduce to granular hydrostatics, 
see Section 2. The hydrostatic problem is fully described 
by three scaled parameters introduced below: the area 
fraction of the particles /, the aspect ratio of the box 
A and the effective hydrodynamic inelastic heat loss pa- 
rameter A 0, By solving the hydrostatic equations 
numerically, we obtain, in Section 2, the spinodal line 
and the critical point of the phase separation in the limit 
of A — ► oo. Section 3 deals with the same limit of A — > oo 
but, in addition, assumes a close proximity to the criti- 
cal point. Here we find the spinodal and binodal (coex- 
istence) lines analytically and determine the structure of 
the steady-state domain wall, separating the more dense 
and less dense stripes in the lateral direction. Effects 
of finite aspect ratio A are addressed in Sec. 4. These 
include suppression of the phase separation by heat con- 
duction in the lateral direction and a change from a su- 
percritical to subcritical bifurcation. Section 5 includes 
a brief discussion of the results and proposes directions 
of future work. 



II. HYDROSTATICS OF PHASE SEPARATION: 
MODEL AND GENERAL RESULTS 

In this Section we formulate the model, briefly review 
the stripe state and compute the spinodal line and the 
critical point of the phase separation for a laterally infi- 
nite system, A — > oo. 



A. Model 

Consider an assembly of inelastically colliding hard 
disks of diameter d and mass m = 1 , moving in a box with 
dimensions L x x L y at zero gravity. Collisions of disks 
with the walls x — and y — ± L y /2 are assumed elas- 
tic. Alternatively, periodic boundary conditions in the 
y-direction can be imposed. We refer to the y-direction 
as the lateral direction of the system. In a laterally in- 
finite system L y — > oo. The particles are driven by a 
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"thermal" wall: a constant granular temperature To is 
prescribed at x — L x . The inelasticity of the particle 
collisions is parametrized by a constant normal restitu- 
tion coefficient r; we will work in the nearly elastic limit 
1 — r 2 <C 1. We also assume a moderate number density 
n: n/n c < 0.5, where n c — 2/{^/id 2 ) is the hexagonal 
close packing density. The last two assumptions allow 
us to employ the Navier-Stokes granular hydrodynamics. 
Throughout this work, we will use for concrete calcula- 
tions the constitutive relations suggested by Jenkins and 
Richman These relations were derived by analogy 

with those obtained in the framework of a successful but 
still empiric Enskog theory (l2|. Detailed comparisons 
with molecular dynamic simulations shows that the er- 
ror margin of the Enskog heat conductivity can reach 
as much as 10 - 15 percent [l^. Still, these relations 
seem to be the best available constitutive relations for 
moderate densities. Importantly, most of the analytical 
results in this work are written in a more general form 
which only assumes a Navier-Stokes structure of hydro- 
dynamic equations. The Jenkins-Richman's relations are 
used only for computing numerical factors. 

Energy input at the thermal wall balances the energy 
loss due to inter-particle collisions, so we assume that 
the system reaches a steady state with a zero mean flow. 
Therefore, the full hydrodynamic equations reduce to hy- 
drostatic equations: 



p = const and V • (kVT) = I , 



(1) 



where p is the granular pressure, T is the granular tem- 
perature, k is the thermal conductivity, and I is the rate 
of energy loss by collisions. Notice that we did not ac- 
count, in the second of Eq. (JJ, for an additional (inelas- 
tic) contribution to the heat flux which is proportional 
to the density gradient. For a dilute gas this term was 
derived in Ref. [l4|- It can be neglected in the nearly 
clastic limit which is assumed throughout this paper. 

The constitutive relations entering Eqs. include the 
equation of state p — p(n, T) and expressions for k and I 
in terms of n and T. In our notation, these relations can 
be written as 0] 
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where v = n(nd 2 /4) is the solid fraction. Let us rescale 
the coordinates hy L x : r/L x — ► r. In the rescaled 
coordinates the box dimensions become 1 x A, where 
A = L y /L x is the aspect ratio of the box. Introducing a 
normalized inverse density z(x, y) — n c /n(x, y) and elim- 
inating the temperature, one can rewrite Eqs. Q as a 



single equation for z(x,y) 0,0, @ : 

V • (F(z)Vz) =AQ(z) 
where F(z) = A(z) B(z), 
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(5) 



which appears in the right hand side of Eq. ©, is the 
hydrodynamic inelastic heat loss parameter. Notice that 
it can be made arbitrary large (by taking large enough 
L x /d), no matter how small the inelasticity q = (1 — r)/2 
is. 

Now we specify the boundary conditions for Eq. 
At the elastic walls x = and y — ±A/2 the normal 
component of the heat flux must vanish. In terms of the 
inverse density z one has V„z = at these three walls. 
Here index n denotes the gradient component normal to 
the wall. Alternatively, for the periodic boundary condi- 
tions we should demand z(x,y + 2ir/A) = z(x,y). The 
constant temperature at the thermal wall x — 1 yields 
the simple condition 



z(x = 1, y) = const 



(6) 



with an a priori unknown constant. As the total number 
of particles N is fixed, the normalization condition 



1 

A 



A/2 dxdy 
-A/2 z(x,y) 



(7) 



should be imposed, where / = (n) /n c is the area fraction 
of the grains and (n) = N/(L x L y ) is the average number 
density of the grains. 

Equation @ with the boundary conditions at the four 
walls and Eq. J7J) make a complete set. Notice that 
the steady-state density distribution is independent of 
the wall temperature To The wall temperature only 
sets the scale of the temperature profile in the system, 
and affects the steady-state pressure. The governing pa- 
rameters of the system are the scaled numbers A, / and 
A. For a laterally infinite system, A — > oo, only two 
governing parameters A and / remain in the hydrostatic 
formulation. 
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Let us define the scaled temperature T 
pressure 

P 



P 



n c Tn 



TU(z). 



T/T and 



(8) 



where 



n(*) = 



l + 2G(z) 



(9) 



Once the steady state density profile and the (uniform) 
steady state pressure p are found, Eq. (jSJ determines the 
steady state temperature profile T[x,y). 

The vector field F(z)Vz entering Eq. J3J is, up to a 
sign, the scaled heat flux. Equation J3J gets simpler if 
we introduce, as a new variable, the scalar potential of 
the heat flux: ip — J z F(z') dz' . To avoid a divergence of 
the integral at infinity, we account for the diverging part 
directly, by extracting the first two terms of expansion of 
F{z) at z — > oo: 



F(z) 
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The integral of f(z) already converges at infinity, and we 
obtain 



1 



iP(z) = -J.Z 



3/2 



+ 
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V*- f(z)dz, (11) 



The potential ip grows monotonically with z, see Fig. ^ 
Now Eq. J2J becomes a nonlinear Poisson equation 



V 2 ip = AQty) 



(12) 



where by Q(ip) we actually mean Q [z(ip)] here and in 
the following. The function Q(ip) is depicted in Fig. [21 
We will be dealing with Eq. (|12f> throughout the paper. 
The boundary conditions for ip{x,y) are identical to the 
boundary conditions for z{x,y): 



0. 



dx 

x—0 

^(1, y) = const , 



(13) 
(14) 



supplemented by either no-flux, or periodic boundary 
condition at the walls y = ±A/2. Note that z is as- 
sumed to be expressed through ip in Eq. (JJJ. 

Importantly, Q(ip) decreases with an increase of ip at 
large enough ip. (In the dilute limit, ip ^> 1, one has 
Q{ip) oc ip^ 1 ^.) This implies non-uniqueness of steady 
state solutions of Eq. I|12f> |l5j|. and opens the way to 
phase separation and coexistence. 



B. Stripe states, spinodal line and critical point 

The simplest steady state of the system is the "stripe" : 
a laterally uniform state which corresponds to the y- 
independent solution z = Z(x) ;2J, or ip — &(x). It is 
described by the equation 



*" = AQ(tf), 



(15) 
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FIG. 1: The effective heat flux potential ip versus the inverse 
scaled density z. The inset shows a blowup of the region of 
1 < z < 4. 
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FIG. 2: The effective heat loss function Q(ip) which appears 
on the right hand side of Eq. Ill 2t . 



with the boundary condition 

tf'(O) = 



and normalization condition 



dx 



o Z[9{x)] 



(16) 



(17) 



Here and in the following the primes denote the x- 
derivatives. As Eq. Ijl5|l does not include the first deriva- 
tive ^'(x), it has "energy" integral and therefore is in- 
tegrablc. A numerical solution however, is more prac- 
tical. Figure [J] gives an example of the density pro- 
file of the stripe state in terms of the scaled density 
n(x)/n c = Z~ x (x) and the auxiliary functions ^(x) and 
Z(x) for A = 344.2 and / = 0.095, corresponding to the 
critical point of the phase separation (see below). 

The stripe state problem I|15 [l -p7 [l can be recast into a 
more convenient initial value problem if we use, instead of 
the normalization condition (|1 7|) . a boundary condition 



*(0) = a. 



(18) 



Indeed, the initial value problem defined by Eqs. i|15[l . 
(|16fl and l|18|) has a unique solution ^(x, a, A). Having 
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FIG. 3: Spatial profiles of the stripe state for A = 344.2 
and / = 0.095, which are the critical values of A and / for 
the constitutive relations (@J. Shown in the main part of the 
figure is the scaled density n(x)/n c = Z" 1 (x). The inset 
shows Z(x) (the solid line) and the heat flux potential 't'(x) 
(the dashed line), obtained by solving Eqs. (1151 - (117| . 



found it, one can calculate the area fraction / = /(a, A) 
from Eq. (|17fl . Importantly, for the Enskog-type con- 
stitutive relations /(a, A) turns out to be a mono- 
tonic function of a for any fixed A. This enables one 
to use the pair of numbers (a, A) instead of (/, A) for 
the parametrization of all possible stripe states. The 
parametrization (a, A) will be often used in the follow- 
ing. 

In a wide region of the parameter space (/, A, A) the 
stripe state undergoes a phase separation instability and 
gives way to a laterally asymmetric state 0, la, S B S 
ll(l | . The instability is driven by negative compressibility 
of the stripe state in the lateral direction [f| 0] , resulting 
from energy loss in particle collisions. Let P = Tli(Z) 
be the scaled steady-state gas pressure of the stripe state. 
In the limit of A — > oo, the spinodal region in the (/, A)- 
plane is defined by the condition (dP/df) A < 0. As P 
is constant in space, it can be conveniently computed at 
the thermal wall x = 1 |6j . Here T = 1 and therefore 

P = P(a,A) = n[zW)}U= n i, a ,A). (19) 
Alternatively, 

1 + 2G[Z X (/,A)] 



P(M) 



Zi(.f,A) 



(20) 



where Z\ = Z(x = 1). The spinodal line of the phase 
separation (again, in the limit of A — *• oo) is determined 
by the condition 



dP(f,A) 
df 



= 0. 



(21) 



Solving Eqs. (fT5|l - l(T^|l . we computed from Eq. f2"U|l the 

P(f) curves at different A, see Fig. These computa- 
tions yield a critical point (/ C ,P C ) [equivalently, (/ c , A c ), 
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FIG. 4: The spinodal line (the solid line) and the P(f) curves 
for the stripe state, for A = 5 ■ 10 3 and A = 280 (the dashed 
and dash-dotted lines, respectively) (a). The two dotted 
straight lines are the dilute-limit asymptotes Pi (A, f) — f 
and P2(A, /) = //2 (a). Figure b shows the P(f) curve for 
A = 5 ■ 10 s in a different scale, to make the region of negative 
compressibility dP/df < more visible. 



or (<z c , A c )]. For A < A c the pressure P increases mono- 
tonically with / (like in an elastic gas), so there is no 
phase separation instability. For A > A c the pressure P 
versus / has a non-monotonic part. The maximum and 
minimum points of P(f) at different A yield the spin- 
odal line. This line in the (/, P)-plane is shown by the 
solid curve in Fig.^Ju Figure 0Ji also shows two straight 
lines. These are the dilute-limit asymptotes of P(f) and 
of the low-density branch of the spinodal line, respec- 
tively (see the next subsection). Figure 0Jd shows, on a 
different scale, the P(/)-dependence for a fixed A within 
the spinodal region, to make the region of negative lat- 
eral compressibility more visible. Figure ^depicts a part 
of the spinodal line, together with the asymptotics of the 
spinodal and binodal (coexistence) lines in the vicinity of 
the critical point, found analytically in Section 3. 

Now let us again use the et-parametrization of the 
stripe solution. Equation (|21fl for the spinodal line is 
equivalent to 



dP(a,A) 
da 



0. 



(22) 
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FIG. 5: The spinodal line (the solid line) and the spinodal 
and binodal asymptotics close to the critical point (the dashed 
lines) in the variables /, P (a) and /, A (b). 



The critical point is the merging point of the maximum 
and minimum points of the P(f) [or P(a)] curve at dif- 
ferent A. Therefore, the critical point is defined by the 
following two conditions: 



dP(a,A) 

da 
d 2 P(a,A) 
da 2 



(o cs A c ) = 0, 



(a c , A c 



0. 



(23) 
(24) 



For an assembly of hard spheres below the freez- 
ing point one has (dp/dn) T > 0, which follows 

(dP/dZ)\ z=Zl (a,A) + Therefore, Eq. (22) for the spin- 
odal line is equivalent to 



a, A) 
da 



0. 



(25) 



while Eqs. (|23|l and l|24|) for the critical point can be 
rewritten as 



(l,o c ,A c ) 



0. 



da 

a 2 * 

— (l,a c ,A c ) = 
da 2 



(26) 
(27) 



(recall that the first argument of function *F and of its 
derivatives stands for x). Solving Eqs. (T5)l. (JTBJl, (fTHfl . 



(|26|l and (|27(l numerically, we find the critical point 



a c 

Ar 



6.580.. 
344.2 . . 



(28) 
(29) 



Using Eqs. (|17() and (|20() . we find the critical point in the 
variables /, P: 



f c = /(o c ,Ac) = 0.0950. 
P c = P(a c ,A c ) = 0.0373. 



(30) 
(31) 



We checked that the maximum density of the stripe 
states, corresponding to the spinodal shown in Fig. 
is less than 0.5 n c , that is within the assumed validity 
domain of the constitutive relations @. At the criti- 
cal point itself the maximum scaled density of the stripe 
state is n(x = 0, a c ,A c )/n c = 1/Z(x = 0,a c ,A c ) = 
0.2086 . . ., a moderate value. 

The critical point, predicted by the Enskog-type gran- 
ular hydrodynamics, agrees fairly well with molecular dy- 
namic simulations by Soto et al. [7J, [10| . The issues of 
accuracy of the hydrodynamic results are discussed in 
Section 5. 



C. Dilute-gas limit 

In the dilute-gas limit, Z > 1, we can obtain the two 
straight-line asymptotes shown in Fig. of the stripe 
pressure P(f, A) and of the low-density branch of the 
spinodal line. Here the stripe state equation is 



dx 2 



Z 3/2 = 3Z~ 1/2 



(32) 



where a different rescaling of the coordinate is intro- 
duced: x = x\J~k. (Notice that in the physical coordi- 
nate Xphys the new rescaling does not include L x .) The 
boundary conditions are 



Z(x = 0) = Z a 



and 



— (x = 0) 

ax 



0. 



where Zq is related to / and A by the normalization con- 
dition J Q Z~ 1 dx — /A 1 ' 2 . The solution to this problem 



Z 



( arccosh\/C + ^C 2 — C 



where C = Z/Z and 
Zo = 



4AV2 



2/AV2 + s inh(2/A 1 /2) 



(33) 



(34) 



In the dimensional units the density profile H33J) is de- 
termined by a single parameter: £ = /A 1 / 2 . The scaled 
pressure of the stripe state is 



P 



1 _ 2£ + sinh(2£) 
~Z[ ~ 4AV2 cosh 2 \ ' 



(35) 
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Z cosh £. Assuming £ <g; 
SJ the dilute-limit asymptote 



where Z\ = Z(x = 1 
1, we obtain from Eq. 

P(f,A) = f- 

The low-density asymptote of the spinodal is deter- 
mined from the condition 

a/A UcA 

This equation, together with Eq. I|35[) . yields coth(£) = £ 
. Substituting this result back into Eq. I)35|l , we obtain 
P(/,A) = //2. 

Using the dilute- limit Eqs. (|33l) and (|34l) . one can ver- 
ify that validity of granular hydrodynamics as an accu- 
rate leading order theory in this problem indeed demands 
nearly elastic collisions (see also Ref. [2()- Equation i|33[) 
implies that the inverse scaled density Z(x) can be pre- 
sented as Z — Zq T{x / Zq) , where function T does not 
include any additional parameters. Therefore, the char- 
acteristic scale of inhomogeneity of the stripe state in 
the a;-direction in the rescaled coordinate x is of the or- 
der of Zq. Going back to the physical coordinate, x = 
A 1 / 2 x p i l ys/ L x ~ q 1 / 2 x p hys/d, we see that the inhomo- 
geneity scale is of the order of Z^d/q 1 ^ 2 ~ (nodq 1 ^ 2 )^ 1 . 
For hydrodynamics to be valid, this quantity should be 
much greater than the characteristic mean free path of 
the particles I which is of the order of (nod)^ 1 . There- 
fore, the validity of hydrodynamics in this system re- 
quires q 1 / 2 <C 1, quite a stringent condition. An addi- 
tional discussion of this condition is presented in Sec. V. 



III. SPINODAL AND BINODAL LINES IN THE 
VICINITY OF THE CRITICAL POINT 

A. Mathematical preliminaries 

In the vicinity of the critical point a c , A c (at A > A c ) 
the two-dimensional, phase separated state is very close 
to the one-dimensional stripe state. Therefore, the spin- 
dal and binodal lines can be obtained by expanding 
il>(x, y) around the stripe-state solution ^{x, a, A) in the 
power series of a — a c and A — A c . For brevity of notation 
the subscripts a and A will denote the partial derivatives 
8/ da and d/dA, respectively, while the prime ' will stand 
for the partial derivative d/dx as before. As will become 
clear shortly, the following functions of x at the critical 
point need to be computed: 



and 



As the stripe solution ^f(x,a,A) is available in quadra- 
ture, the derivatives of \I/ with respect to parameters a 
and A are known. More practical, however, is a differ- 
ent approach. As shown in Appendix, all these functions 
[and an additional function $(x) that we will need, see 
Eq. HtjUfl below] can be expressed as solutions of the linear 
differential equation 



with different source terms S(x) and different bound- 
ary conditions. For the first derivatives, and ^> a , the 
source term vanishes. Therefore, the rest of the functions 
can be expressed through and ^o, see Appendix. 

B. Spinodal line in the vicinity of the critical point 

It follows from Eqs. lf2U|) and (|2"7|) that the first non- 
vanishing term in the expansion of ^(x, a, A) in the pow- 
ers of a— a c at the critical point is the cubic term (a— a c ) 3 . 
Therefore, we should keep the following terms in the ex- 
pansion: 

a, A) = *(x, a c , A c ) + A c V A (x, a c , A c ) 5 



+a c $ a (x, a c , A c )u + y ^ aa (x, a c , A c ) u 2 



+ a c A c ^ a A(x, a c , A c ) uS + -f- ^ aaa {x, a c , A c ) u 3 , (37) 

6 

where an order parameter u — a/a c — 1 and control pa- 
rameter S = A/A c — 1 have been introduced. 

To obtain the equation of the spinodal line we differ- 
entiate Eq. (|37|1 with respect to a, put x = 1 and use 
Eqs. l|25 Jl -ip7 )l . The result is 



or 



where 



, a 2 c ^ aaa (l,a c , A c ) 2 

o = ; ; ; — r U 

2A c * aA (l,a c , A c ) 

A - A c = Ax (a - a c ) 2 , 

*aaa(l,a c , A c ) 



A x = - 



2* aA (l,a c ,A c ) ' 



(38) 



(39) 



(40) 



One can see from Eq. (|38|l that, close to the critical 
point, 6 = 0(u 2 ). That is why we could neglect the term 
proportional to * A A 5 2 = 0{u A ) in Eq. (|37jl . 

The coefficients 'i'aaai^, a. c , A c ) and \E' a A(l, a c , A c ) can 
be computed numerically, see Appendix: 



*„ a „(l,a c ,A c ) = 0.02434... , 
* aA (l,a c ,A c ) = -0.001926.. 



(41) 
(42) 



so Ax = 6.317. . .. 

Eq. H39f) can be rewritten in terms of A and /. Using 
Eq. I|17|) . we expand Z(x 1 a, A) near the critical stripe 
solution Z c (x) = Z(x, a c , A c ) up to the first order in a—a c 
(which suffices close to the critical point) and obtain 



f-fc = fa(a c , A c ) (a - a c ) 



where 



w"(x) - A c [9(x, a c , A c )] w(x) = S(x) 



(36) 



/a(a c ,A c ) = - / - 

J \Z c (x 



* a (x,a c , A c ) 



J [Z c (x)Y F(Z c (x)) 



■ dx . 



(43) 



(44) 
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Evaluating this integral numerically, we obtain 
f a (a c , A c ) = -0.004396 .... As the result, Eq. ^ can 
be rewritten as 



A - A c = A 2 (/ - f c f , (45) 



where 



A 2 = A 1 [f a {a c ,A c )] 2 = 3.268...-10 5 . (46) 

Now we compute the spinodal line in the /, P-plane. 
Expanding Eq. I|19l) in the vicinity of z = Z(l, a c , A c ), 
we obtain 



j P = P c + nf ) * A (l,a c ,A c )(A-A c ) + ... , 
where the higher-order terms are negligible, and 

t(c) _ 



(47) 



n 





dTL\ 


{f(z) 


dz J 



= -2.969...- 10" 



z=z[*(l,a c ,A c )] 



(48) 

(c) 

The negative value of 11^ ' is a consequence of our defi- 
nitions of z and ip and of the condition (dp/dn)T > 0. 
Combining Eqs. (|47fl and l|45|) we obtain: 



P-Pr 



-A 5 (f~f c ) 2 



(49) 



where A$ = — ^"^(l, a c , A c )n^. A numerical calcu- 
lation (see Appendix) gives ^(l, »c, A c ) = 0.2203..., 

therefore A5 = 21.38 The spinodal asymptotics l)49|) 

and (|45l) arc shown, together with the full spinodal line, 
in Figs. [5]} and[5p, respectively. 



C. Two-phase coexistence and binodal line in the 
vicinity of the critical point 

1. Laterally non-uniform states 

For a steady state with a broken lateral symmetry the 
function ip(x, y) depends on its two arguments. In a later- 
ally infinite system the asymptotics of ip( x j 2/) at 7/ — > ±00 
correspond to two different stripe states. Therefore, it is 
natural to replace the no-flux or periodic boundary con- 
ditions in the lateral direction by the condition 



ip(x,y^> ±00) = ip±(x) . 



(50) 



where ip-(x) ^ tj)+(x). One way to solve the prob- 
lem (I12ll - (|14|) and (|50|l is to introduce an unknown func- 
tion a(y) so that 



ip(0,y) = a(y) , 
a{y — ► ±00) = a± = const 



(51) 



with a_ 7^ a + . What equation should a(y) satisfy close 
to the critical point? Here we can look for ip(x, y) in the 
form of a weakly and slowly modulated stripe state: 



ip(x, y) = ^[x, a(y),A] + <fi(x, y) . 



(52) 



where a(y) = a c [l + u(y)] is a slow function of y, <fi(x, y) 
is a small correction to "J, and a y-dependent order pa- 
rameter u(y) <C 1 is introduced. We will see shortly that 
the characteristic length scale of a(y) (the domain wall 
width) is of the order of <5 -1 / 2 ~ it -1 , while (f> ~ u 3 . 
Hence, every y-derivative introduces smallness of order 
u. Making the Ansatz (|52|l in Eq. I|12|) and neglecting 
terms of a higher order than it 3 , we arrive at the follow- 
ing linear problem for (f>(x,y): 



dl(j>- A c Q^ c {x)} (f> = -a c V a (x, a c , A c ) , (53) 



<K0,y) = £W(0,y)=0. (54) 

The first boundary condition in Eq. 154|) clmiminates the 
arbitrariness in the choice of 4>{x, y), while the second one 
follows from Eqs. I|13|l and (|16fl . Additional boundary 
conditions include 

a, A) + 0(1, y) = const (55) 
at the thermal wall [see Eq. 1|14|) ]. and 

4>(x,y — > ±00) = (56) 

[see Eq. Q5U[l]. Notice, however, that Eqs. (|55|) and 
(|56|l do not enter the problem (|53|) and (|54|l for cj)(x,y). 
Equation (|53|l can be solved by separation of variables: 



(x,y) = a c ${x) 



d 2 u(y) 
dy 2 



(57) 



The function $(ir) is the solution of the following prob- 
lem: 

$" - A c Q^[* c (x)] $ = -* (<r, a c , A c ) , (58) 

$(0) = $'(0) = 0, (59) 

which again belongs to the class of equations (|36|l . The 
solution (see Appendix) is 



<f>(x) 



A c Q(a c 



X 

^ a (x,a c ,A c ) J y'^ a d£ 



X 

*'(x,a c ,A c ) J *ld£ 



(60) 



where the functions 'J'' and ^> a under the integrals over £ 
have arguments £, a c and A c . Now we impose the bound- 
ary condition (|55|l : 

cPu 
] dy~ 2 



*(l,a,A) + ^$(1)3-^ = const, (61) 
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As the second term in the right side of Eq. I|55|) is of order 
u 3 , we should keep terms up to u 3 in the expansion of 
iff (x, a, A) near the critical point a Cl A c . This expansion 
has the same form as Eq. (|37[1 . with the only difference 
that now u depends on y. Evaluating this expansion at 
x = 1 and using the definitions of the critical point [Eqs. 
(|26[) and i|27[l ]. we arrive at the desired equation for u(y): 



\ S — A3 u 3 + A4 —— - = a = const , 
dy z 



where 

A 3 = - 



a^ ann (l,a c , A c ) _ ajAi 
6A c *oA(l,a c ,A c ) ~ 3A C 



= 0.26485 . 



A 4 = 



*(1) 



A c * aA (l,a c ,A c ) 
and [see Eq. JSJJ] 



= 0.5212. 



$(!) = - 



V(l,a c ,A c) , , 



A c Q(a c ) 



*f dx = -0.3457 . 



(62) 



(63) 



(64) 



(65) 



2. Domain wall and the binodal line 

One integration of Eq. (|62|l yields: 
2 a / r \ 2 



A A 



M 
2 



u 2 - 



+ 2au + /3 . 



(66) 



where (3 = const. To find the constants a and /3, we are 
using the boundary conditions in the lateral direction. In 
an infinite system these are 



u(y — > ±00) = u± = const 



(67) 



with w_ 7^ u_|_. Equations (|66|l and (|67|l appear in numer- 
ous problems of domain wall structure, see for example 
Ref. |l6j. The boundary conditions l|67|l yield a = (3 = 0, 
so the domain wall solution is 




tanh 



2A 4 



(y - yo) 



(68) 



where ± refers to two possible orientations, and the ar- 
bitrary constant yo describes the position of the domain 
wall. The domain wall solution exists only if 5 > 0, 
that is A > A c . Eq. I|68|l confirms our assumption that 
the characteristic width of the domain wall is C(<5 -1 / 2 ). 
Returning for a moment to the physical (unsealed) vari- 
ables, we see that the domain wall width is 0(L X / '5 1 / 2 ) 
which is much larger than L x . 

The values of the order parameter far from the domain 
wall arc 



u± 



a± 

Cic 



1 



(69) 



Equation (|69|l . rewritten in terms of a and A, defines the 
binodal (coexistence) line: 



A - A c = ^- (a - a c ) 2 , 



or, in terms of / and A, 



A- A,, 



(70) 



(71) 



Compare these expressions with Eqs. (|39H and (|45l) for 
the spinodal. 

The binodal line can be also expressed in terms of / 
and P. The derivation almost coincides with that for 
the spinodal line, see the end of Sec. IIIIBI The only 
difference is that now we substitute into Eq. I|47|l the 
binodal relation (|71() rather than the spinodal relation 
|@SJ. The result is 



P-Pr 



A, 



U-fcf 



(72) 



The physical meaning of the binodal asymptotics l|72() 
and H71JI is straightforward. Firstly, the two coexisting 
stripes with f — fi and / = / 2 have equal pressures: 
= P{h)- Secondly, close to the critical point, fi 
and /2 are symmetric with respect to f c , that is {f\ + 
/ 2 )/2 = f c - The binodal asymptotics (f?^ and (fTT|) are 
depicted in Figs. and Eb, respectively. 

Note that our results for the spinodal and binodal 
lines [see Eqs. lj4*5*|) and (J7J)] are consistent with the 
results obtained by Soto et al. 0, El in the framework 
of their phenomenological "van der Waals equation" . In- 
deed, our Eqs. (|43|l and (|7T)> have the same quadratic 
forms as those obtainable from the van der Waals equa- 
tion @>EJ- Furthermore, the ratio of the coefficients of 
the spinodal and binodal lines l|45|l and l|71|) equals three, 
in agreement with what the van der Waals equation pre- 
dicts. (Note that the coefficients themselves of the van 
der Waals equations have not been derived yet.) In ad- 
dition, Soto et al. 01 reported a binodal line found in 
molecular dynamics simulations of this system for a mod- 
erate value of the parameter e — d/ L x = 0.01. See Sec. 
5 for a discussion of the role of this parameter. 

Soto et al. also suggested an interpretation of the 
binodal line in terms of the "Maxwell's construction". 
By analogy with the classical van der Waals gas [l7j . 
Maxwell's construction can be written as 



P(/,A)-P(/!,A) 
P 



df = 0, 



(73) 



where the factor f 2 in the denominator comes from equal- 
ity dV = d(l/f) = -df/f 2 , where V = 1// is the 
(scaled) specific volume. How does Eq. (|73|) compare 
to our result (|72|l ? Consider expression l|19fl for P(a, A). 
Close to the critical point, it can be expanded in the 
powers of u and 5: 



1 

P = P C + A C P A 5 + a c A c P aA uS + -a 3 c P ( 
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, (74) 
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where all the derivatives of P are evaluated at the critical 
point. Going over from u to / — f c , one obtains 

Ja 

+ ^(/-/c) 3 + -.., (75) 
6/1 

where the numerical coefficient f a is given by Eq. (|44|l . 
and the dots stand for terms of the order of (/ — f c ) A 
and higher. At i5 > (when phase separation occurs), 
the coefficient in front of / — f c is negative, while that in 
front of (/— f c ) 3 is positive. Note also that, in the vicinity 
of the critical point, S = 0(f — f c ) 2 ■ Now we substitute 
Eq. I|75|) . with the higher-order terms neglected, into Eq. 
(|73|l . In this order of perturbation theory one should 
put / = f c in the denominator of the integrand. As 
the result, Eq. I|73[l reduces to the simple relation (/i + 
/2)/2 = f c obtained above. Of course, the reason for 
this agreement is the closeness to the critical point. In 
this sense, the Maxwell's construction (|73|) does not give 
anything new. Moreover, many different constructions, 
for example, 

h 

J[P(f,A)-P(f 1 ,A)}df = 0, (76) 

fx 

are equally applicable in the vicinity of the critical point. 
Of course, the Maxwell's construction would be valuable 
if it were shown to be true far from the critical point. 
At present there is no reason to believe this is the case, 
and the theoretical form of the binodal line far from the 
critical point remains unknown. 



IV. FINITE SIZE EFFECTS: MARGINAL 
STABILITY AND BIFURCATIONS 

A. Marginal stability surface 

The results presented in Sections 2 and 3 are valid in 
the limit of A — > 00. Already in the first paper on 
the phase separation instability it was found that the in- 
stability is suppressed when the aspect ratio of the box 
A = L y /L x is less than a threshold value A*(/, A). The 
physical mechanism of suppression is heat conduction in 
the lateral direction which tends to erase the lateral tem- 
perature (and, therefore, density) inhomogeneity. How 
to generalize the spinodal line, obtained for A — > 00, to 
finite A-s? Let us consider a small sinusoidal density per- 
turbation, in the lateral direction, around a stripe state. 
The fastest growing (or the slowest decaying) perturba- 
tion is the one with the longest wavelength, compati- 
ble with the boundary conditions in the lateral direction 
0- ■ Consider a three-dimensional parameter space 
(/, A, A) and define in it a two-dimensional marginal sta- 
bility surface !F{f, A, A) = 0. By definition, at any point 



on this surface the growth rate of the longest pertur- 
bation is equal to zero. The marginal stability surface 
represents a natural generalization of the spinodal line. 
Importantly, this definition reduces to that of the spin- 
odal line as A — > 00 [see Ref. 0] and Eq. (|8"TI) below]. 
The marginal stability surface can be computed by lin- 
earizing Eq. (|12fl around the stripe state and solving the 
resulting linear eigenvalue problem. Calculations of this 
kind were done previously for large A far from the crit- 
ical point H HI. Three typical cross-sections of the 
marginal stability surface are shown in Fig. As ex- 
pected, the instability region shrinks as A goes down. 
As the result, the critical point moves toward smaller f-s 
and larger A-s as A decreases. The monotonic depen- 
dence of the critical point position on A is explained by 
the monotonic increase of the lateral heat conduction as 
A decreases. 

The threshold value of the aspect ratio A* (A, /) has a 
minimum at some / 0,0. The respective minimum value 
A m i„ depends only on A. In this work we performed a 
systematic investigation of this dependence. We found 
that A m i n goes down monotonically as the parameter 
A — A c is positive and increases. Though this monotonic 
decrease looks like a power law in the log- log plot, see 
Fig. [7^, it is actually not. Figure [7|d shows the same 
dependence on a different scale. Two different asymp- 
totes are clearly seen. The first of them, at A > 1, 
was obtained previously: A m i„(A) = AA -1 / 2 , where 
A = 52.14... 0, 0. In this regime the eigenfunction 
of the marginal stability problem is exponentially local- 
ized at the elastic wall x = 0. The second asymptote is 
valid close to the critical point A = A c , where A m i„ di- 
verges: A min = 42.085 ... (A - Ac) -1 / 2 . This asymptote 
is derived analytically in the next subsection. 

An additional interesting issue is the change of bifurca- 
tion character, predicted by a weakly nonlinear analysis 
of the steady state problem close to the marginal stabil- 
ity surface. At fixed A, the bifurcation is supercritical 
on an interval /-(A) < / < /+(A) which lies within the 
spinodal interval (/i,/?)- On the intervals f\ < f < f_ 
and f+ < f < fi the bifurcation is subcritical 0, 0]. 
The next subsection addresses the finite size effects in 
the vicinity of the critical point, where everything can be 
calculated analytically. 



B. Steady states and bifurcation types close to the 
critical point 

This subsection addresses two-dimensional steady 
states in a laterally finite system close to the critical 
point. As A min diverges at the critical point [see Eq. i|97)|) 
below], we assume that A, though finite, is very large. 
The starting point here is the same Eq. I|62|) . but on 
a finite interval —A/2 < y < A/2, so we replace the 
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FIG. 6: Cross-sections of the marginal stability surface by 
three planes: A = oo (the solid line), A = 1.788 (the dash- 
dotted line) and A = 0.523 (the dotted line). The dashed lines 
show the dilute-limit asymptote /A 1 ' 2 = 1.1997..., obtained 
analytically [fj, and the spinodal and binodal asymptotics 
near the critical point. 



boundary conditions i|tj7|) by the no-flux conditions 

= . (77) 



du 
dy 



du 
dy 



Periodic boundary conditions can be treated in a similar 
manner. Integrating Eq. (|62|) over y from —A/2 to A/2 
and using the boundary conditions l|77|l . we determine a 
and rewrite Eq. (|1)2"|) as 



(78) 



Ai — + u 5 - A 3 u 3 - (u S - A 3 u 3 ) = . 
dy z x ' 

where (...) denotes spatial averaging: 

A/2 

(...) = ! J (...)dy. 

-A/2 



Introduce the rescaled coordinate ?/* = y/A, order pa- 
rameter u, = uA and control parameter (5* = A 2 6. 
Equation l|78|) keeps its form in the rescaled variables, 



A, 



dru* 



A 3 ul-(u*)6* + A 3 (ul)=0, (79) 



while the boundary conditions become 



du* 



, > du* , 
■1/2) = —(,,. 



1/2) = 0. (80) 



As the aspect ratio A drops out in these variables, a 
universal description can be obtained. 

Obviously, any y-independent state u* = const solves 
the problem l|79|l - H80() (a y-independent state is nothing 
but a stripe state). What is the condition for the ap- 
pearance of (weakly) y-dependent solutions? When a 
y-dependent solution does appear, what is the type of 
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FIG. 7: The dependence of the threshold value of the as- 
pect ratio A for the suppression of the phase separation on 
the parameter A — A c (a). The phase separation instabil- 
ity occurs when A > A m j n (A). Figure b shows that the 
apparent "straight line" in figure a actually includes two dif- 
ferent asymptotes and a crossover between them. The two 
asymptotes are shown by the dashed lines: the A > 1 asymp- 
tote, A m i„(A) = 52.14... A _1,/2 , and the asymptote near 
the critical point, A m i„ = 7r (A c A4) 1 / 2 (A — A c )~ 1//2 , where 
A4 — 0.5213 . . .. The results of numerical marginal stability 
analysis are denoted by circles. 



bifurcation? To address these questions, we seek for a 
weakly y-dependent solution in the form 



it* = (u*) + di sin 7ry* + a 2 cos 2ny* 



(81) 



We substitute Eq. ljST)l into Eq. (|7Tj|l and treat the terms 
originating from a 2 cos 27ry* as small corrections. Ex- 
panding up to 0(a 3 ), we obtain the following two alge- 
braic equations: 

3 

5. - n 2 A 4 - 3A 3 (u*) 2 - - A 3 a\ + 3A 3 (u„) a 2 = , 



- A 3 (u*) a\ + (<5» - 47r 2 A 4 - 3A 3 {u»)) a 2 = . (82) 

Putting here a\ = a 2 = 0, we obtain the marginal stabil- 
ity condition 



3A 3 ul + 7T 2 A A , 



(83) 
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where we have omitted the spatial averaging of u% , as 
it becomes trivial on the marginal stability curve. The 
marginal stability curve is shown, as the thick solid line, 
in Fig. [S] Now we consider nonzero amplitudes a\ and 
02 in Eq. I|82l) and eliminate a 2 in favor of a\ . Above the 
marginal stability curve (|83|1 . but close to it, we obtain 
the following equation for the bifurcation curve: 



(84) 



3A 3 4 V~ tt 2 ^4 

One can see that, on the marginal stability curve, the 
bifurcation is either supercritical (when the term in the 
parentheses in the right hand side of this equation is pos- 
itive), or subcritical (when this term is negative). The 
change of character of the bifurcation occurs at the points 

1/2 

/, : ) = ±| I =±3.116... , (85) 



2A, 



ir 2 A 4 



12.86. 



(86) 



which lie on the marginal stability curve Q83[) . see Fig. [8] 




FIG. 8: The marginal stability curve for a laterally finite sys- 
tem (the thick solid line) is shown together with the spinodal 
and binodal lines for an infinite system (the thin solid line 
and the thin dotted line, respectively). This is made possible 
by using the rescaled order parameter u* and control param- 
eter <5» . The two black circles show the points 1851 and 1861 , 
where the bifurcation changes its character from the super- 
critical (between the circles) to the subcritical (outside). 

Note that the spinodal and binodal lines for an infinite 
system can be expressed through and <5* : 5* = 3 A 3 u 2 
and (5* = A 3 u% , respectively. This makes it possible to 
show all the curves (the spinodal and binodal for an in- 
finite system and the marginal stability curve for a finite 
system) on the same plot, see Fig. EI 

Going back from the rescaled variables (5* and m* to A 
and a, we obtain for the marginal stability curve: 



A - A c = Ai (a - a c ) 2 + 
or, in the variables /, A: 

A - A c = A 2 (/ - f c ) 2 + 



7T 2 A 4 A C 

A 2 

7T 2 A 4 A C 



(87) 



(88) 



The physical meaning of this result becomes transparent 
when one compares it with Eq. 145(1 for the spinodal in 
the infinite system. 

Solving Eq. (|88|) for A, we can find the threshold value 
of the aspect ratio A* at given / and A: 



A,„ 



[A - A c - A 2 {f - fc) 2 ] 1 ^ 2 



(89) 



The phase separation instability occurs at A > A*. One 
can see that A* diverges on the spinodal line. The min- 
imum value of A* at a given A corresponds to / = f c : 



A _ kaTkI" 



(A-A C )V2 



(90) 



This quantity diverges at the critical point, as the lower 
asymptote in Fig. shows. 

Finally, the two points of change of the bifurcation 
character, [Eqs. I|85|l and (|86|l ]. form a line, 



A -A, 



■A 1 (a- a c y 



in the plane a, A, which corresponds to the line 



A-A c = -A 2 (/-/ c ) 2 , 



(91) 



(92) 



in the plane /, A. This line lies within the spinodal in- 
terval of the problem corresponding to A — > 00. 



V. SUMMARY AND DISCUSSION 

In this work we have employed granular hydrostat- 
ics to determine the phase diagram of a prototypical 
driven granular gas which exhibits spontaneous symme- 
try breaking and van der Waals-like phase separation. 
We determined the spinodal line and the critical point 
of the phase separation. We computed the spinodal and 
binodal (coexistence) lines close to the critical point . Ef- 
fects of finite lateral size of the confining box have been 
also addressed. These include determining the line of 
change of the bifurcation character from supercritical to 
subcritical. 

The shape of the binodal (coexistence) line far from 
the critical point is still unknown. To handle this prob- 
lem within the framework of the hydrostatic theory one 
needs to solve the nonlinear Poisson equation l|12|) in a 
laterally infinite box. Most likely, this can only be done 
numerically, in a sufficiently long box. It is clear, how- 
ever, that the binodal line in the variables (a, A) is solely 
determined by the function Q(ip) which encapsulates all 
the necessary information about the equation of state, 
heat conductivity and inelastic energy loss. 

It is worth mentioning that our qualitative results (the 
phase separation instability, the existence of the criti- 
cal point and spinodal and binodal lines, the suppres- 
sion of the instability by the lateral heat conduction and 
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the change of bifurcation character) are robust and in- 
dependent of the small details of the constitutive rela- 
tions. Furthermore, most of our analytical results are 
presented in quite a general form which only assumes a 
Navier-Stokes structure of the hydrodynamic equations 
(in the limit of nearly elastic collisions). We used the 
Enskog-type constitutive relations [111 ] only for comput- 
ing numerical factors (and for the computations far from 
the critical point). Obviously, the accuracy of our quan- 
titative results cannot be better then the accuracy of the 
Enskog-type constitutive relations. Therefore, a 10-15 
percent error margin should not be surprising. 

There is an important additional aspect that can be- 
come crucial when comparing the hydrodynamic theory 
with molecular dynamic simulations. As was already 
mentioned, for the hydrodynamics to be quantitatively 
accurate, the mean free path of the particles must be 
much less than any hydrodynamic length scale. It is well 
known that, for a system with a thermal wall, the leading 
correction to hydrodynamics enters the boundary condi- 
tion in the form of a temperature jump, proportional to 
the ratio of the mean free path and the characteristic 
hydrodynamic length scale 0, ^| . Indeed, within the 
Knudsen layer, whose size is comparable to the mean free 
path, hydrodynamics is inapplicable, while the particle 
velocity distribution significantly deviates from Maxwell 
distribution: the temperature of the incoming particles 
is less than that of the outgoing particles [as the outgo- 
ing particles have their normal velocity randomized ac- 
cording to a Maxwell distribution with a (fixed) higher 
temperature]. The effective temperature at the wall, for 
the purpose of a hydrodynamic description in the bulk 
(that is, outside the Knudsen layer), is always less than 
To 0,03, as is indeed observed in MD simulation of this 
system 2j. As the result, the pressure of the system is 
reduced (notice that the density does not change in this 
order; it can change only in the next, Burnett order). 
The pressure reduction will obviously cause a shift of the 
critical point. For our numerical results to be sufficiently 
accurate the parameter e = d/L x must be very small, 
so that the mean free path is indeed much less than the 
system dimensions. As A ~ qs~ 2 = const (for example, 
at the critical point A = A c ), a very small e implies an 
extremely small inelasticity q. Despite the severe limita- 
tion intrinsic to it, the nearly elastic case is conceptually 
important, just because hydrodynamics is supposed to 
give here a quantitatively accurate leading order theory. 

An interesting direction for future work is the phase 
separation dynamics. Here the hydrostatic equations 
should give way to the full set of hydrodynamic equa- 
tions. Close to the critical point, however, a reduced 
description of the dynamics should be possible. It has 
been suggested that such a description is provided by the 
"van der Waals equation" 0, LUj. Though the van der 
Waals equation does capture qualitatively much of the 
phenomenology of the phase separation as seen in the 
MD simulations 0, EJ; its systematic derivation from 
the equations of granular hydrodynamics is still lacking, 



and the coefficients of the normal form are yet unknown. 
Importantly, our hydrostatic results close to the critical 
point, encapsulated in Eq. fully agree with those 

predicted from the van der Waals normal equation. Fur- 
thermore, Eq. (|62|l provides quantitative relations be- 
tween the yet unknown coefficients of the van der Waals 
equation. 

Another avenue of future work requires going beyond 
hydrodynamic description, as it includes two types of 
fluctuations in this system. The first of them was ob- 
served in molecular dynamic simulations inside the spin- 
odal region at A = 11050 (that is, far from the critical 
point) in a wide region of aspect ratios around the finite- 
size threshold value A* 0. It was found that these fluc- 
tuations dominate the dynamics of the system, so they 
were called giant fluctuations. The second type of fluc- 
tuations is expected to occur, by analogy with the classi- 
cal van der Waals phase transition, in a close vicinity of 
the critical point. The fluctuations in these two regimes 
should be describable in the framework of "fluctuating 
hydrodynamics" of Landau and Lifshitz 19], general- 
ized to granular gases in the limit of nearly elastic col- 
lisions. Fluctuating hydrodynamics is a Langevin-type 
theory which takes into account the discrete character 
of particles by adding delta-correlated noise terms in the 
momentum and energy equations |l9j| . The fluctuations 
appear in this approach as a hydrodynamic response of 
the system to the Langevin noise. Unfortunately, this 
exciting direction of work is hindered by the fact that 
the Langevin term, which accounts for the discreteness 
of the inelastic energy loss in the energy equation, has 
yet to be calculated Q. 

Finally, it was assumed throughout the paper that the 
granulate is driven by a thermal wall. In experiment a 
rapidly vibrating wall is usually used. Though qualita- 
tively similar, the phase diagram of the case of a rapidly 
vibrating wall can be quantitatively different 0, |H • 
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APPENDIX. COMPUTING THE COEFFICIENTS 

Consider the stripe solution ^(x, a, A). The expan- 
sions in the vicinity of the critical point, used through- 
out the paper, include several derivatives of this function 
which need to be evaluated. These are ^ aa , 
^Aj ^Aa, ^aaa, and additional function 4> [see Eqs. 
and (|59|l ]. all of them evaluated at the critical point 
a = a c , A = A c . One can easily show that each of these 
functions is a solution of the linear problem 

w"(x) - A c (tf c ) w(x) = S(x) (93) 
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w(x) 


S(x) 


Initial conditions at 


x = 






w(0) 


iw'(O) 


f'(x,a c ,A c ) 








Ac Q (a c ) 


■tya(x,a c , Ac) 





1 





* A (a;,a c ,A c ) 


Q(* c ) 








$oa(l,a c ,A c ) 


Ac Qv>y> (* c ) #a 








^/Aa(x,a c , Ac) 


Ac Q^V (*c) *A *a + Qi, (*c) *a 








^aaa(x,a c ,A c ) 


3A C QvV> (* c ) *a*aa + A c (*c) *a 








$(x,a c ,A c ) 











TABLE I: The source terms S(x) and the initial conditions at x — for Eq. 19311 . for each of the auxiliary functions shown in 
the first column. The functions in the first and second columns have the same arguments. 



with different source terms S(x) and different initial con- 
ditions at x = 0. The source terms and initial conditions 
are listed in Table. [I] 

Let us show, as an example, the derivation of S(x) 
and of the initial conditions for two of the functions: \l/ a 
and ^Aa- The starting point is Eq. (|15fl with the initial 
conditions ((Tfj|l and (|T%)) for the stripe solution ^(x, a, A). 
The stripe solution depends on a and A because they 
enter either the equation, or the initial conditions. Let 
us differentiate the both sides of Eq. I|15fl with respect to 
a. We obtain 



<'(z,a,A)-A^ 



i/i=>P(x,a,A) 



* a (x,a,A) =0. (94) 



Once ^(x, a, A) is known, Eq. (|94(l is a linear differential 
equation for ty a . Now let us differentiate, with respect 
to a, Eqs. I|16(l and (|18fl . We obtain the relations 



* tt (0, a, A) = 1 and %(0, a, A) = . 



(95) 



which serve as the initial conditions for the same function 
\& . As we see, the source term S(x) vanishes in this case. 

Differentiating Eqs. I|94|) and l|95(l with respect to A, 
we arrive at the following problem for the function ^ a A- 



< A (z,a,A)-A 



dip 



ip=<l>(x,a,A) 



^ aA (x,a,A) 



dQ 

dip 



i/j=*(x,a,A) 



t> a (x,a,A)+ 



< A (0,a,A) = 0. 



(98) 



Again, once \l/ a and \I/a are known, Eq. l|9rj|l is a linear 
equation for ^ a \. Here the source term is nonzero, while 
the initial conditions (|97fl and i|98|) are zero. Substituting 
in Eqs. ll9"4l - (|9"Sl) a — a c and A = A c , we obtain the 
second and fourth rows of Table [Q The other rows of 
Table 1 can be obtained in a similar way. 

The functions ^'{x, a c , A c ) and ^ a {x,a c , A c ) are spe- 
cial, as they satisfy the homogeneous form of Eq. I|93|l : 
S(x) = 0. Therefore, once they are found, the rest of the 
functions from Table 1 can be expressed through them: 



w(x) — 



1 



A c Q(a c ) 



*\x) J * o (0S(0*- 





C 1 ^ a (x) + C 2 ^'(x), (99) 



where G\ and C% are integration constants. To satisfy 
the zero initial conditions at £ = 0, we must choose 
C\ = C% = in all cases. Furthermore, evaluating w(x) 
at the thermal wall x = 1, we observe that the term pro- 
portional to ^a(x) in Eq. (|9"9"|l vanishes, as ^(1) = at 
the critical point. Therefore, for all functions from Table 
1, except iff' and we obtain 



A 



d 2 Q 

dip 2 



ip='$(x,a,A) 



V A (x,a,A.)V a (x,a,A), (96) 



w(l,a c ,A c 



A c Q(a c 



(100) 



* aA (0,a,A) =0, 



(97) 
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